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AN ACCURATE STRAPDOWN DIRECTION COSINE ALGORITHM 


By John W. Jordan 
Electronics Research Center 


SUMMARY 

This report presents an algorithm for computing the direction 
cosine matrix of a strapdown inertial system* It is designed for 
efficient computation on a general purpose digital computer. It 
includes a correction for commutativity error and has no trunca- 
tion error. Its high accuracy is demonstrated by simulation. 

The simulation technique can be applied to the evaluation of other 
direction cosine algorithms as well, permitting a comparative 
evaluation . 


INTRODUCTION 

Strapdown inertial systems involve two computational pro- 
cesses which are not present in conventional inertial systems 
employing a gimbal structure. One is associated with the trans- 
formation of vehicle acceleration from body coordinates to navi- 
gational coordinates : 


a N (t) = C T (t) a B (t) (1) 

where 

a D = vector acceleration of the vehicle in body coordinates 
B 

a N = vector acceleration of the vehicle in navigational 
coordinates 

C = the 3x3 matrix of direction cosines which represent the 
transformation from body to navigational coordinates. 

The second process is the determination of the direction 
cosine matrix by the real-time integration of the matrix differ- 
ential equation: 


C(t) = Ji(t)C(t) 


( 2 ) 



where ft(t) is a skew-symmetric matrix of vehicle rates measured 
in body coordinates : 


ft (t) = 


' 0 
-w z (t) 
W y (t) 


w z (t) 

0 

-w (t) 

x 


-w y (t) 

w x (t) 

0 


Both Eqs , (1) and (2) can be solved by a general-purpose computer 

they are replaced by finite difference equations. The frequen- 
cy with which the equations must be solved is determined by a 
system accuracy specification and the anticipated environment; 
that is, by the vehicle acceleration and angular rate profile. 
Generally, the solution of Eq. (1) is not too demanding of the 
available computational resources (ref. 1). Eq. (2), on the 
other hand, may require a substantial percentage of the computer's 
time because of the inherent complexity of the numerical methods 
used and the high iteration rates which may be required to achieve 
the specified accuracy. This report is concerned with the solu- 
tion of Eq. (2) on a general-purpose digital computer. In general, 
the airborne computer will calculate a matrix B(t) which only 
approximates the true cosine matrix C(t). 


FINITE DIFFERENCE EQUATION 

The computer will employ a finite difference equation for 
propagating the cosine matrix: 

B* = $B (3) 

where 


B* is the new solution at step n 
B is the old solution at step n-1. 

A suitable transition matrix is given in reference 2: 


$ = I + + k 2 e 2 


where 0 is the skew-symmetric rotation matrix: 
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0 d> 
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(4) 


(5) 
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Let: 


2 2 2 2 

c h =<b +<p +<b . 

o x y z 

The scalar parameters k^ and k 3 are given in Table I along 
with a scalar k 3 which will be used subsequently. 

TABLE I 


r — 1 

k 2 

k 3 

sin d> 

o 

1 - cos d> 
Y o 

tan (<|> o /2) 

o 

-e- 

cO 2 

Y o 

o 

-e- 


Equation (4) has appeared frequently in the literature 
(refs. 1, 4, 5, 7) as an analytical solution of the differential 
Eq. (2) for a special form of the vehicle rate matrix ft(t). For 
that particular case, the rotation matrix (Eq. (5)) is simply 
the integral of fi(t). In this report, Eq. (4) is used in a 
different way than in the previous literature. It is presented 
as a general, exact expression for describing the rotation of a 
cosine matrix. The problem then becomes the determination of an 
appropriate rate matrix (Eq. (5)), so that the difference equation 
(Eqs. (3) and (4)) represents the solution of the direction cosine 
problem. In general, the rotation matrix is no longer equal to 
the integral of 0, {t) , although it must, of course, reduce to this 
form for the special case when Eq. (4) does represent the solution 
of Eq. (2) . 

Before considering the determination of the required rotation 
matrix, Eqs. (3) and (4) will be put in a form suitable for effi- 
cient computation on a general=purpose computer. 


COMPUTATIONAL ALGORITHM 


An efficient computational algorithm may be ‘devised by 
expressing the matricies in terms of vectors. The direction 
cosine matrix requires three vectors. 


B = 




b 


3 


3 




The rotation matrix, being skew- symmetric , can be represented by 
one vector. 


4> = 


The computational algorithm is as follows: 


<f> Q = 


(6) 


Compute 


for 


k l' k 3 
= k 1 4> 


*3 = V 

j = 1/2 

a . = b . x d>T 
3 3 1 

3. = a. x cK 


bj = bj + (a j + 3..) 

* * * 
b 3 = bi x b 2 


The last step is based upon the fact that the direction 
cosine matrix is orthogonal. The calculation of and k^ is 
discussed in an Appendix. The following section will discuss 
the determination of the rotation vector 4> . 


DETERMINING THE ROTATION VECTOR 

Strapdown inertial systems measure vehicle rotation by 
gyroscopes mounted ("strapped down") to the frame of the vehicle. 
Reference 3 provides an expression for the rotation vector in 
terms of the integral of the body rates sensed by the gyroscopes 

t 

<p =f (j) ( X ) dX + a (7) 

o 


4 


where 


cp = The rotation vector 

w = The vehicle rate vector sensed by the system gyroscopes 

a = A vector equal to the area traced out by each body axis 
on a unit sphere. 

Most inertial systems use rate-integrating gyroscopes which 
provide the integral of the vehicle rates: 

h 

A =f cu (A ) dA 
o 

The quantity A is the vector output of the system gyroscopes 
over a sampling interval h. 

The area vector, a, is found (ref. 3) by applying the two 
dimensional Green's theorem: 

ir h 

a = ~ 2 J 4> x <f> dA 
o 

Then Eq. (7) becomes 

l 

A = <j ) + -jJ<j>x < t>dA (8) 

o 

The problem now is to determine <j> from the gyroscope output 
vector A. Figure 1 illustrates the idea of a pre-processor. 

The gyroscopes provide the vector. A, which is sampled at a 
frequency f s . The pre-processor calculates <f> from A and uses 4> 
to update the direction cosines at a frequency f c . The pre- 
processor algorithm has two constraints - it must provide suffi- 
cient accuracy and it must be computationally simple. 


PRE-PROCESSOR ALGORITHM 

One possible pre-processor algorithm can be obtained from 
Eq. (8) by considering the cross-product term to be small com- 
pared to the other right-hand term. This is equivalent to choos- 
ing the sampling frequency, f s , high enough so that the direction 
of the rotation vector changes little over the sampling period. 
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This assumption results in the approximation 


<f> = A - ijj- 1 Axu> dA (9) 

*x> 

In order to evaluate the integral, it is necessary to obtain 
the vehicle rate vector, uj, from its integral. A, which is pro- 
vided by the system gyroscopes. This can be done by fitting a 
polynomial to a sequence of A vectors and then analytically dif- 
ferentiating the polynomial to obtain the vector ui . Reference 3 
introduced the concept of a pre-processor and a polynomial approx- 
imation in order to obtain an estimate of the vehicle angular 
rate matrix, ft(t), given a sequence of gyroscope outputs. The 
angular rate matrix was then used to integrate Eq. (2) directly, 
using a Runge-Kutta algorithm. In this report, the pre-processor 
and polynomial approximation are used to determine the equivalent 
Euler angle rotation from the gyroscope outputs. Reference 6 
discusses the use of a pre-processor to provide increased accuracy 
for direction cosines calculated by the method of quaternions. 
However, the pre-processor algorithm is equivalent to Eq. (11) of 
this note; the more general form of Eq. (10) is not derived or 
simulated. Reference 5 also derives a correction term equal to 
Eq. (11) by a completely different method. 

Let 


A (A) = d Q + d ± \ + d 2 X 2 

with (f) = dg being the accumulated value of the rotation vector 
since the last calculation of the direction cosines. The two 
remaining coefficients can be determined by two gyroscope samples 
A^ and See Figure 2. The pre-processor algorithm is then 

determined by substution to be 

ip = A + A 2 (10) 

* 12 
<J> = <p + ip + -^<pxip + j- (A^ x A 2 ) 

The vector, <j> , is updated every two gyro samples or at a frequency 
f s /2. This high-speed calculation is continued until the direction 
cosines are updated at a lower frequency, f c . The vector <f> is 
then set equal to zero and the high-speed calculations are re- 
sumed. For the case where f c = f s /2, the initial value of <f> will 
always be zero and the pre-processor algonithm is simply 

<fp = A^ + A 2 + (A^ x A 2 ) (11) 
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COSINE ERRORS 


For some applications, a simplified form of Eq. (4) is 
sufficient. Approximations for and k 3 can be used instead of 
the true values of Table I. The resultant error may be called a 
truncation error. The truncation error is completely under the 
control of the designer and it will vanish when the calculations 
of k! and k 3 are made exact. For the remainder of this note, 
it will be assumed that the accurate values of Table I are used 
and that there is no truncation error involved in updating the 
direction cosines. However, the Appendix will consider trunca- 
tion error resulting from the approximate calculation of k-^ and k 3 . 

There is another error source which results in what is often 
called commutativity error. For the algorithm presented in this 
note, it is caused by the approximate nature of the pre-processor 
algorithm. In general, the pre-processor provides an estimate 
$ of the required rotation vector <p . Although the commutativity 
error could be defined as the difference between the vectors <j> 
and $, in this note it will be defined as the erjror in the 
direction cosines which results from the use of <f> instead of the 
truerotation vector (J> . See Figure 3. One special case is 
important. If 


<P x <j> *= 0 


then Eq. (8) reduces to 


4> - A 

and the pre-processor is error free. This condition can also be 
expressed as 


a) x A = 0 (12) 

which is satisfied if the vehicle rate vector, w, has a fixed 
axis in space over the sampling interval. If this condition is 
satisfied, there is no commutativity error and the direction 
cosines can be determined exactly. When Eq. (12) is satisfied, 
the vehicle rate matrix, Q(t), corresponding to the vector co is 
said to be "self-commutative" (ref. 5). For this special case, 
the rotation vector is simply the gyroscope output vector. 


SIMULATION 

The effectiveness of the pre-processor can be demonstrated 
by simulation. In order to do this, it is necessary to have an 
analytical solution for the direction cosines with a non— self- 
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commutative rate matrix. Although no general solution is known, 
a class of test solutions are easily constructed as follows. Let 


C 


c c 


where 


C = fiC 


C A * VA 


C B = Vb 


If fi A and fig are self-commutative, then C A , Cg and hence C 
are easily found. The angular rate vector for C is easily shown 
to be 


W 


“a + C A "b 


and 


-t + h 

A = A. + J C a) dX 
A J A B 

Next, u A and oj b are chosen so that the integration can be done 
analytically. A simple choice of 



1 

• 

cr> 

J 


1 

o 

“a = 

l 

0 o 

1 

W B = 

o a\ 


was used. This gives a vehicle rate vector of 


U) 


.6 

.6 cos.6t 
- . 6 sin . 6t . 


a severe coning motion. The commutativity error of Figure 3 can 
be expressed as a single scalar parameter (ref. (4). 

1 

e = k f trace (E T E )1 2 
c L c c J 
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where k = (57.3) (60)/t 

s 

t = simulation time (min) 
s 

This error parameter has the units of degrees per hour. A value 
of t s = 20 minutes was used for all simulations. An additional 
parameter N was defined by 


f = N f 
s c 

so that N is the number of gyro samples before the direction 
cosines are updated. The results for N = 1 correspond to a 
direct application of Eq. (6) every time the gyros are sampled 
without using a pre-processor. Since the accurate values of 
and k 3 of Table I were used, none of the algorithms have any 
truncation error. The results are plotted as a function of f s and 
f c in Figure 4 and 5. Since updating the direction cosines in- 
volves the bulk of the computation, the plot of the commutativity 
error vs. f c most closely approximates a plot of error vs. com- 
puter loading. All algorithms except N = 1 involve the calcula- 
tion of the pre-processor every two gyro samples. For larger 
values of N, this requires proportionally more computation and 
this fact should be kept in mind when interpreting Figure 5. The 
solution with N = 2 has an additional advantage since its pre- 
processor algorithm (Eq. (11)) is less complex than that of the 
others (Eq. (10)). 


CONCLUSION 

The pre-processor significantly reduces the secular commu- 
tativity error with the algorithm N = 2 providing the most accu- 
racy for an equivalent computer loading. Eqs . (6) and (11) pro- 

vide a high-accuracy algorithm for general-purpose computers which 
has low commutativity error and no truncation error. 

Since the truncation error is easily analysed (Appendix A) 
and is completely under the control of the designer, it follows 
that the commutativity error represents a more significant test 
of an algorithm's preformance. Quite possibly, a more effective 
pre-processor can be designed, perhaps by using a higher order 
polynomial approximation. However, it is not the inherent 
accuracy of an algorithm which is important but rather its reali- 
tive effectiveness at a specified level of computer loading. 
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APPENDIX A 


TRUNCATION ERROR 

Truncation error results from approximations made in calcu- 
lating k-L and k 3 . This error may also be determined by simulation 
as in Figure 3. For a self-commutative rate matrix, ref. (5) 
provides an analytical solution for the truncation error. If the 
approximations for k^ and k 3 are of low accuracy, the resulting 
cosine matrix may become non-orthogonal and an orthogonality 
correction can be introduced. The correction, however, is com- 
plex and the added computational time would be more advantageously 
used in improving the accuracy of k^ and k 3 . Reference 5 con- 
tains an analytical solution for the truncation error when an 
orthogonality correction is used. It also considers the error 
introduced by finite computer word length. Refs. (1) and (7) also 
contain an analysis of these error sources. 


CALCULATION OF kj^ k 3 

The calculation of kj^ and k 3 may be done using standard 

proceedures for trigonometric functions. However, the expres- 
sions in Table I cannot be used since they may require division 
by zero. For some applications a very simple approximation may 
be sufficient. For example, refs (1), (5), (6), and (7) discuss 

kf = 1, k 2 = 0, 1/2, and 1/3. However, by using the pre-processor 
to correct the commutativity error and by calculating more accu- 
rate values of k^ and k 3 , all the calculations can be done at a 
lower frequency. An excellent discussion of the calculation of 
trigonometric functions is contained in ref. (8) . What follows 
is based on that reference. 


The calculation of tan (x/2) provides a good basis for all 
the trig functions. For example: 


sin x 


tan tx/2) 

1 + tan 2 (x/2 1 


Hence 


k l = 


2 2 
1 + *0 V 
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and it remains to calculate k^. Let P n (x) the polynomial 

n 

p n (x) = a 3 x3 
y=o 

The tangent function can be calculated as 

2 

tan (ttx) = X P n (X ) 

with the coefficients for a nth order approximation given in 
ref. (8). In order to determine k 3 , only P n (x 2 ) is necessary. 
This will be a polynomial in 4 >q 2, which is provided in Eq. (6). 
The polynomial may be calculated by Horner's nested form 


v k ■ x 2 v k + i + a k k ■ "- 1 ' n ~ 2 ' "• 0 

p n (x 2 ) - v o 

which requires n additions and n multiplications. Reference 8 
discusses a streamlined form as well, which requires only 3 
multiplications and 7 additions for a 6th order polynomial. 

(The streamlined form was not used in the simulations reported 
here.) The calculation of P n (x 2 ) may also be done in a sub- 
routine which is used for all trigonometric calculations in the 
airborne computer. 


13 



REFERENCES 


1. Wilcox, J. C.: A New Algorithm for Strapped-Down Inertial 
Navigation. IEEE Transactions on Aerospace and Electronic 
Systems, vol . AES-3, no. 5, September 1967. 

2. Crandall, S. H. et al: Dynamics of Mechanical and 
Electromechanical Systems. McGraw-Hill Book Company, New 
York, 1968. 

3. United Aircraft Corporate Systems Center: A Study of the 
Critical Computational Problems Associated with Strapdown 
Inertial Navigational Systems. NASA CR-968, April 1968. 

4. Sullivan, J. J. : Evaluation of the Computational Errors of 
Strapdown Navigational Systems. AIAA Journal, vol. 6, 

no. 2, February 1968. 

5. Jordan, J. W.: Direction Cosine Computational Error. NASA 
TR-R-30 4 , March 1969 . 

6. TRW Systems: Body-Fixed, Three-Axis Reference System Study, 
Phase II, Final Report. Prepared for NASA George C. 

Marshall Space Flight Center, Huntsville, Alabama, under 
Contract No. NAS8-20209, 15 December 1966. 

7. Morgan, A. M. : Computational Errors in the Generation of 
the Direction Cosine Matrix in a Strapdown System. Journal 
of the Institute of Navigation, vol. 14, no. 1, Spring, 1967. 

8. Hart, J. F. et al. : Computer Approximations. John Wiley 
and Sons, New York, 1968. 


14 


NASA- Langley, 1969 21 C”80 


